REPORT  DOCUMENTATION  PAGE 

- — - AFRL-SR-BL-TR-01-  - 

Public  reporting  burden  for  this  collection  of  information  is  estnated  to  average  1  hour  per  response,  including  the  time  tor  reviewing  instruct  mg  and  reviewing 

the  collection  of  information.  Send  comments  regarding  this  burden  estvnat^  or  any  other  aspect  of  this  collection  of  information,  mclud  ^  ^  7  J  tor  Information 

Operations  and  Reports,  1215  Jefferson  Oavis  Highway,  Suite  1204.  Arlington,  VA  22202-4302.  and  to  the  Office  of  Management  and  Budg 


1.  AGENCY  USE  ONLY //Mb'S 

2.  REPORT  DATE 

31  Oct  00 

3.  F  _ _ 

Final  Tech  Report  1  Mar  96  to  28  Feb  98 

4.  TITLE  AND  SUBTITLE 

Theoretical  and  Numerical  Validation  of  Scalar  EM  Propagation  Modeling  Using 
Parabolic  Equations  and  the  Pade  Rational  Operator  Approximation 

5.  FUNDING  NUMBERS 

F49620-96- 1-0060 

6.  AUTHOR(S) 

Dr.  Ronald  Brent 

7.  PERFORMING  ORGANIZATION  NAME{S)  AND  ADDRESS(ES) 

Department  of  Mathematical  Sciences 

University  of  Massachusetts  Lowell 

Lowell,  MA  01854 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research 

AFOSR/NM 

801  N.  Randolph  Street,  Rm  732 

Arlington,  VA  22203-1977 

lO.SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

F49620-96- 1-0060 

11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  unlimited. 

13.  ABSTRACT  (Maximum  200  words} 

The  problem  associated  with  the  absorbing  layers  had  been  addressed  and  solved.  The  code  has  been  benchmarked  against 
analytical  solutions  and  shown  to  be  accurate  to  within  .01  dB.  The  enclosed  figures  show  excellent  agreement  between  the 
computed  and  exact  solution  over  a  variety  of  frequencies.  While  the  transmission  loss  for  an  infinite  medium  is  frequency 
independent,  which  is  why  all  the  solutions  are  identical,  the  construct  of  the  absorbing  layer  is  frequency  dependent.  This  is 
why  I  checked  a  variety  of  input  frequencies  between  1  MHz  and  50  MHz.  The  following  pages  contain  output  showing 
exact(anaiytical)  solutions  and  2  numerical  solutions,  one  with  the  absorbing  layer  and  one  without.  One  can  clearly  see  that 
the  absorbing  layer  correctly  dampens  any  energy  associated  with  reflections  from  nonphysical  (finite  numerical  domain) 
boundaries. 


Unfortunately  it  took  most  of  the  contract  period  to  do  so,  hence  none  of  the  original  objectives  were  accomplished. 


Also  appended  to  this  report  is  an  updated  user  manual  containing  user  information,  as  well  as  an  interpolation  and  source 
sensitivity  analysis. 

D’HC  A,  4 


14.  SUBJECT  TERMS  “  20010124  137 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 

OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


15.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


20.  LIMITATION  OF 
ABSTRACT 

_ UL 

Standard  Form  2S8  (Rev.  2*89)  (EG) 

Prescribed  by  ANSI  Sid.  239. 1 6 

Designed  using  Perform  Pro.  WHS/DIOR,  Oct  94 


UNCLASSIFIED 


UNCLASSIFIED 


UNCLASSIFIED 


Theoretical  and  numerical  validation  of  scalar  EM  propagation  modeling  using 
parabolic  equations  and  the  Fade  rational  operator  approximation. 

F49620-96- 1-0060 


Final  Report 


Dr.  Ronald  Brent 

Department  of  Mathematical  Sciences 
University  of  Massachusetts  Lowell 
Lowell  MA  01854 
(978)  934-2440  (Phone) 
(978)934-3053  (Fax) 

Ronald_Brent(^uml.edu  (email  after  6/1/2000) 


ORIGINAL  OBJECTIVES 


*  Benchmarking  the  numerical  techniques  used  including  a  sensitivity  analysis. 

*  Developing  and  enhancing  the  numerical  code  by; 

*  Studying  how  to  interpolate  the  refraction  profiles  on  the  field  mesh. 

*  Implementing  a  finite  element  algorithm. 

*  Incorporating  troposcatter  and  sky  wave  propagation 

INTERMEDIATE  OBJECTIVES 

*  Build  a  working  absorbing  layer  below  the  Earth's  surface. 

*  Benchmark  the  layer  against  simple  analytical  models 

FINAL  STATUS 

The  problem  associated  with  the  absorbing  layers  had  been  addressed  and  solved.  The 
code  has  been  benchmarked  against  analytical  solutions  and  shown  to  be  accurate  to  within  .01 
dB.  The  enclosed  figures  show  excellent  agreement  between  the  computed  and  exact  solution 
over  a  variety  of  frequencies.  While  the  transmission  loss  for  an  infinite  medium  is  frequency 
independent,  which  is  why  all  the  solutions  are  identical,  the  construct  of  the  absorbing  layer  is 
frequency  dependent.  This  is  why  I  checked  a  variety  of  input  frequencies  between  1  MHz  and 
50  MHz.  The  following  pages  contain  output  showing  exact(analytical)  solutions  and  2 
numerical  solutions,  one  with  the  absorbing  layer  and  one  without.  One  can  clearly  see  that  the 
absorbing  layer  correcth  dampens  any  energ\’  associated  with  reflections  from  nonphysical 
(finite  numerical  domain)  boundaries. 

Unfortunately  it  took  most  of  the  contract  period  to  do  so.  hence  none  of  the  original 
objectives  were  accomplished. 

.'\lso  appended  to  this  report  is  an  updated  user  manual  containing  user  information,  as  well  as  an 
interpolation  and  source  sensitivity  analysis.  ' 


Absorbing  Layer  Comparison  (f=lMHz,  hmaxl==3000=-hmax2)  wf=15 


Numerical  Solutiorr  (with  absorbing  layer) 


Absorbing  Layer  Comparison  (f=5MHz,  hmaxl=-hmax2=3000)  wf=2 


(ap)  “TO 


Numerical  Solution  (with  absorbing  layers 


Absorbing  Layer  Comparison  (f=10MHz,  hmaxl=-hmax2=3000)  wf=l 


(ap)  ureo 


Numerical  Solution  (with  absorbing  layers 


Absorbing  Layer  Comparison  (f=20MHz,  hmaxl=3000=-hmax2)  wf=l/3 


Numerical  Solution  (with  absorbing  layers) 


Absorbing  Layer  Comparison  (f=25MHz,  hmaxl=3000=-hmax2)  wf=l/3 


(OP)  UIBO 


Range  (km) 


Range  (km) 


SSPiUser  Manual 


by 


Dr.  Ronald  Brent 

Department  of  Mathematical  Sciences 
University  of  Massachusetts  Lowell 
Lowell  MA  01854 
(978)  934-2440  (Phone) 
(978)934-3053  (Fax) 
Ronald_Brent@uml.edu 


Abstract 


This  report  contains  updated  user  information  for  the  program  SSP .  This  program  is  an 
implementation  of  the  Fade  series  approximation  for  parabolic  equations  modeling  one-way 
scalar  2-D  electromagnetic  propagation.  The  code  can  be  run  on  PC’s  with  Pentium  processors. 
The  current  program  is  set  for  20,000  notch  points  in  height  and  NPADE  <  8 .  The  code  is 
written  in  generic  FORTRAN  and  has  been  compiled  using  Microsoft-FORTRAN  Powerstation. 
The  accompanying  CD  contains  source  and  executable  codes  for  SSP  as  well  as  EPADE,  the 
code  responsible  for  calculating  the  Fade  coefficients.  Also  on  the  CD  is  a  graphics  program 
GRAPH  1  and  SSPPREP  an  input  file  preparation  program  dsigned  to  make  setting  up  input 
files  more  user  friendly.  GRAPH  is  a  screen  graphic  contour  programs  for  displaying  results 
immediately.  This  program  is  written  in  C  and  source  codes  are  also  included. 


1.  Program  Description 

Under  certain  simplifying  assumptions  one  may  reduce  the  solution  of  the  full  3-D 
vector  electromagnetic  propagation  problem  to  that  of  two  simpler  scalar  models.  These 
problems  correspond  to  horizontally  or  vertically  polarized  wave  fields.  A  full  discussion 
of  the  analysis  and  derivations  leading  to  such  equations  is  given  in  References  1  and  2. 
The  primary  assumption  leading  to  this  reduction  or  splitting  of  the  vector  problem  into 
scalar  pieces  is  medium  symmetry  in  at  least  one  coordinate  direction.  The  reader  is 
referred  to  Ref  1  before  reading  this  report.  We  will  assume  that  the  user  is  familiar  with 
that  report  and  the  notations  and  descriptions  contained  therein.  The  program  SSP  is  the 
numerical  implementation  of  the  results  of  that  report  It  is  a  fimte  difference  approach  for 
solving  parabolic  equations  modeling  scalar  EM  propagation. 

The  program  SSP  consists  of  a  main  routine  which  calls  a  total  of  13  other 
subroutines  or  functions.  The  subroutine  hierarchy  is  illustrated  in  Figure  1 .  Subroutine 
Descriptions  are  as  follows: 

INPUT  This  subroutine  is  responsible  for  reading  in  all  input  parameters  other  than 
field  profiles.  This  input  resides  in  the  first  7  lines  of  the  file  SSP.IN,  and 
the  next  n  lines  determining  the  terrain  elevation.  If  the  user  requests  an 
absorbing  layer,  INPUT  calls  the  subroutine  SPONGE,  and  uses  the 
function  CRVCORR. 

CRVCORR  This  function  is  used  in  implementing  the  curvature  correction  when  using 
the  Earth-Flattening  transformation. 

PROFL  Responsible  for  input  of  all  profile  data  which  resides  in  the  lines  following 
the  topography  input.  PROFL  calls  the  subroutine  READQ. 

READQ  Responsible  for  reading  in  medium  profiles  and  interpolating  the  profiles 
over  the  entire  height  grid. 

CC  Responsible  for  computing  the  width  of  the  artificial  absorbing  layer  at  the 

top  of  the  atmosphere.  For  a  given  firequency,  using  tabular  data,  absorbing 
layer  widths  are  interpolated  and  passed  back  to  INPUT. 
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GREEN  Creates  a  Green’s  function  starting  field.  It  uses  the  function 
FGS(ARG)  =  exp(  -ARG  /  4) 

MODES  Creates  a  Homogeneous  Normal  Mode  Starting  Field.  This  is  a  nice  startup 
field  because  at  allows  limitation  of  the  aperture  using  THMAX. 

GAUSS  Creates  the  typical  Gaussian  Starting  field.  It  uses  the  function 
FGR(ARG)  =  (1 .4467  -  0.4201  *ARG)*exp(  -ARG  /  3.0512) 

MATRC  This  subroutine  sets  up  the  matrices  used  in  solving  the  linear  system 
associated  with  the  Fade  discretization  of  the  Parabolic  Equation. 

UPDAT  This  subroutine  is  used  to  update  the  matrices  whenever  there  is  a  change  in 
the  bottom  topography. 

SOLVE  This  subroutine  solves  the  tridiagonal  system  of  equations  associated  with 
the  discretization.  This  subroutine  is  called  every  time  a  step  in  range  is 
taken. 

The  program  flow  is  quite  simple.  For  a  given  medium,  the  basic  flow  is  a  loop  in 
range  in  which  the  solution  of  a  tridiagonal  matrix  problem  is  calculated.  After  each  step, 
receiver  loss  values  are  interpolated  and  then  written  to  a  file.  The  matrices  involved  in  the 
linear  system  are  composed  of  tridiagonal  entries  dependent  upon  grid  width,  dz,  medium 
parameters,  such  as  refraction  and  conductivity,  as  well  as  interface  coupling.  Each  time 
any  one  of  these  parameters  changes  the  matrices  must  be  re-computed.  The  two 
subroutines  MATRC  and  UPDAT  are  responsible  for  setting  up  these  matrices.  The 
subroutine  MATRC  is  used  whenever  there  is  a  new  set  of  profiles  in  the  input  run-stream. 
For  terrain  elevation  changes  UPDAT  is  used.  The  subroutines  UPDAT  and  MATRC  do 
exactly  the  same  thing  however  when  there  is  only  a  change  in  terrain  elevation,  and  not 
refiuction,  one  need  only  re-compute  the  matrix  elements  for  a  small  sub-set  of  equations 
in  the  entire  linear  system.  The  program  runs  far  more  efficiently  using  UPDAT.  Running 
the  code  using  MATRC  alone  is  possible,  and  in  fact  is  a  good  check  if  adapting  this  code 
in  that  area. 

A  general  flow  of  the  program  is  given  in  Figure  2.  This  is  only  a  rough  outline  of 
the  code  to  give  the  reader  a  feel  for  the  program. 
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SSP  PROGRAM  STRUCTURE 


Figure  1 
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SSP  Flow  Chart 


Figure  2 
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cont. 


Figure  2  cont. 
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11.  Input  Run-Stream 


The  input  for  the  program  comes  from  two  files.  The  first  file  EPADE.DAT  must  be 
produced  by  using  the  program  EPADE.  This  program  automatically  reads  physical 
information  from  the  primary  input  file,  SSP.IN,  and  produces  EPADE.DAT  containing 
the  Pade  coefficients.  The  input  file  SSP.IN  is  where  the  user  specifies  all  the  information 
pertaining  to  souTce  and  receiver  configurations  as  well  as  the  medium  topography  and 
refraction  profiles.  The  disk  containing  the  source  and  executable  codes  also  has  a 
directory  with  some  sample  input  files.  It  is  recommended  that  a  user  first  become  familiar 
with  those  files  and  modify  them  for  their  own  use. 

The  input  in  the  file  SSP.IN  is  grouped  into  three  sections.  The  first  section,  which  is 
contained  in  the  first  seven  lines  of  the  file,  is  where  the  source  and  receiver  configuration 
is  specified.  Output  options  are  also  specified  there  along  with  parameters  needed  for  the 
screen  graphics  post-processing  program  GRAPH.  The  next  section  is  the  terrain  profile. 
The  remaining  section  itself  consists  of  sub-sections.  The  medium  dielectric  and 
conductive  properties  are  entered  as  a  series  of  range-independent  sets  of  profiles.  That  is 
at  specified  ranges,  height  dependent  profiles  are  entered.  There  is  no  interpolation  in 
range.  The  medium  is  assumed  vertically  stratified  until  the  following  range  with  changes 
is  reached.  A  complete  description  of  the  input  records  v^th  associated  definitions  is  given 
in  the  following  pages.  The  next  section  contains  some  example  input  run-streams  and 
their  associated  output. 
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Input  Run-Stream 

************ 

CARDl 

************ 

TITLE 

TITLE;  Title  for  Input  Data  and  Graphics  Plot.  After  the  Title  (20  Characters) 
one  should  type  the  word  "TITLE",  i.e.  in  columns  21-25.  This  is 
necessary  for  the  post-processing  graphics  program,  GRAPH. 

************ 

CARD  2 
************ 

FREQ,  ZS,  ZR,  POLFLG,  ABSRBFLG 

FREQ:  Source  Frequency  in  MHz. 

ZS:  Source  Height  (above  MSL)  in  meters. 

ZR:  Receiver  Height  (above  MSL)  in  meters. 

POLFLG:  Polarization  Flag  (  POLFLG=0  for  Horizontal  Polarization  POLFLG=l 

for  Vertical  Polarization  ) 

ABSRBFLG:  Atmospheric  Absorbing  layer  Fiag  (  ABSRBFLG=1  will  extend  the 

atmosphere  by  the  amount  ^'^DTH  as  computed  by  the  subroutine  CC.) 

************ 

CARDS 

************ 

RMAX,  DR,  NDR 

RMAX:  Maximum  Range  for  Calculations  in  kilometers. 

DR:  Computational  range  step  in  meters. 

NDR;  Number  of  notch  point  skips  in  range  in  writing  entire  field. 
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************ 


CARD  4 

************ 

HMAXl,  HMAX2,  DZ,  NDZ,  ZTPLT,  ZBPLT 

HMAXl :  Maximum  height  of  Atmosphere  in  meters. 

HMAX2:  Maximum  Depth  of  Terrain  in  meters. 

DZ:  Computational  height  step  in  meters. 

NDZ:  Number  of  notch  point  skips  in  height  in  writing  entire  field. 

ZTPLT:  Used  for  contour  plot.  Top  of  the  graph  in  meters 

ZBPLT:  Used  for  contour  plot.  Bottom  of  the  graph  in  meters 

NOTE:  All  heights  are  measured  relative  to  MSL  (z=0).HMAX2  gives  the 

elevation  of  the  bottom  of  the  calculational  domain.  (Usually  negative) 

************ 

CARDS 

************ 

GL(i)  i=l,13 

GL(I):  These  are  contour  level  values  (positive)  for  use  by  the  post  processm^ 

program  GRAPH,  which  produces  screen-graphic  contour  plots  of  loss 
values 

************ 

CARD  6 

************ 

NPADE,  (IFLAG(I),I=1,4) 

NP  ADE:  The  number  of  terms  in  the  Pade  approximation  series. 

IFLAG(l):  Write  Flag  for  the  output  file  PROFL.OUT.  This  file  contains 

interpolated  input  profiles.  IFLAG(l)  =  1  turns  the  write  statements  on. 


IFLAG(2):  Write  Flag  for  the  output  file  TRNPRF.OUT  which  contains 

interpolated  terrain  profile  heights.  IFLAG(2)  =  1  turns  the  wnte 
statements  on. 

IFLAG(3):-  Write  Flag  for  the  output  file  SSPl  .OUT.  IFLAG(3)  =  1  turns  the  write 
statements  on.  If  using  the  graphics  program  GRAPH  this  must  be  set 
equal  to  1 . 

IFLAG(4):  Output  option  flag  for  artificially  setting  loss  values  =  500  dB  for 

receiver  locations  below  the  terrain.  This  will  make  for  more  distinct 
graphics  output  when  using  the  post-proccessing  program  GRAPH. 
IFLAG(4)=1  causes  500  dB  substitution.  If  choosing  this  option,  it  is 
best  to  pick  GL(13)  =  499. 

CARD? 

************ 

ISTRT,  THMAX 

ISTRT;  Starting  field  type  flag 

=  I  gives  a  Gaussian  starting  field. 

=  2  gives  a  Green's  fimction  staring  field. 

=  3  gives  a  Normal  Mode  starting  field  using  a  point  omni-directional 
point  source.  Aperture  is  limited  by  THMAX 

THMAX:  Maximum  angle  for  mode  starting  field.  Useful  in  limiting  the  source 

aperture.  Used  only  in  ISTRT=3. 
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************** 


CARD  8+  Terrain  Elevation  Data  Input  ( <  500  data  points) 
************** 


RD(1),  D(l) 
RD(2),D(2) 


Range  RD  in  kilometers  and  elevation  D  in  meters. 


RD(N1),  D(N1) 


*************** 

CARD  8+Nl  Atmospheric  N  Profile. 
*************** 

Z(l),  CW(1) 

Z(2),  CW(2) 

Height  Z  in  meters. 

Refiraction  CW  in  N-units. 

NOTE:  Z(i)  must  be  less  than  Z(i+1). 
Z(N2),  CW(N2) 


******************* 

CARD  8+N1+N2  Terrain  Dielectric  Input  Profile. 
******************* 

Z(l),  CB(1) 

Z(2),CB(2) 

Height  Z  in  meters. 

Relative  Permittivity  CB. 

NOTE:  Z(i)<Z(i+l). 

Z(N3),  CB(N3) 


*****:fc** ************** 


CARD  8+N1+N2+N3  Terrain  Conductivity  Input  Profile. 
********************** 

Z(l),COND(l) 

Z(2),  COND(2) 

Height  Z  in  meters 
Conductivity  COND  in  MHO/m 
NOTE:  Z(i)  must  be  less  than  Z(i+1). 

Z(N4),  COND(N4) 


************************** 

CARD  7+N1+N2+N3+N4 
************************** 


RP 


RP :  Range  in  kilometers  of  the  next  profile. 

iliJtc************  ************ 

CARDS  7+N1+N2+N3+N4-H- 

Repeat  the  previous  three  sections  for  each  range  profile. 

Once  an  input  file  has  been  set  up,  the  user  then  runs  the  program  EPADE.  This 
program  automatically  reads  the  input  data  DR,  FREQ,  and  NPADE,  from  the  file  SSP .IN 
and  then  produces  a  file  called  EPADE.DAT  which  contains  the  complex  series 
coefficients.  This  file  is  then  read  by  the  program  SSP.  The  progrzun  EPADE  is  a 
modification  of  one  written  by  Dr.  Michael  Collins.  A  description  of  the  numerical 
methods  used  in  that  program  can  be  found  in  Ref.  3.  You  will  get  a  message  fi'om 
EPADE  if  the  code  does  not  converge.  If  this  is  the  case  you  must  either  decrease  the 
range  step-size  DR,  or  increase  the  number  of  terms  in  the  series,  NPADE. 


Program  SSPREP 


1  have  written  a  precursor  program  called  SSPPREP  which  is  designed  to  create 
input  files  for  the  novice  user.  It  prompts  the  user  for  all  physical  mformation,  such  as 
source  and  receiver  configuration,  terrain  elevation,  and  atmospheric  makeup,  and  then 
creates  an  input  file  called  SSP.IN  in  the  format  appropriate  for  the  code  SSP. 

The  code  SSPPREP  allows  for  terrain  following  refractive  profiles.  I  put  this  in 
since  there  is  some  data  available  that  suggests  that  near  coastal  regions,  atmospheric 
refraction  profiles  may  tend  to  slope  up  with  the  land  as  the  terrain  rises.  Terrain 
topography  can  be  input  as  (x^)  data  pairs  (typical),  or  as  a  sinusoidal  model  with  user 
specified  amplitude  and  period.  Both  the  input  files  used  in  samples  1  and  3  were  created 
using  SSPPREP. 

If  you  use  SSPPREP  the  input  file  that  is  created  should  have  parameter  values  that 
should  yield  convergent  solutions.  There  is  a  possibility  of  the  code  crashing  if  the 
parameters  ZTPLT  and  ZBPLT  are  not  chosen  correctly. 

III.  Output  Options 


There  are  several  output  options  for  the  code.  Primary  output  is  in  the  two  files 
SSP  LOUT  and  SSP2.0UT.  For  a  requested  receiver  height,  SSP2.0UT  contains  the  loss 
values  at  each  ra~ge  step  between  r  =  0  and  r  =  UMAX.  The  file  SSPl.OUT  is  a  binary 
file  used  by  the  program  GRAPH.  It  contains  field  loss  values  on  a  grid  specified  by  the 
user  in  Card  4  through  parameters  ZTPLT  and  ZBPLT.  These  are  the  elevations  of  the  top 
and  bottom,  respectively,  of  the  graph.  NDR  and  NDZ  give  the  output  steps.  Currently  the 
graphics  programs  accept  grid  sizes  with  the  total  number  of  range  points  between  250  and 
1000,  and  the  total  number  of  height  points  between  125  and  500.  If  the  user  requests  too 
many  height  output  points  the  program  will  automatically  change  NDZ  to  accommodate 
the  output.  Range  steps  NDR  must  be  correctly  entered  by  the  user  however,  if  the 
graphics  program  is  used.  Choose  NDR  so  that  RMAX/(NDR’'‘DR)  is  between  250  and 
1000,  where  UMAX  is  measured  in  meters.  1  wished  the  output  to  this  file  to  be  as  general 
as  possible  for  other  uses.  This  is  why  NDR  and  NDZ  were  not  hard  coded. 
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For  debugging  purposes,  it  was  also  nice  to  have  the  interpolated  medium  refraction 
values  for  the  entire  height  array.  By  choosing  IFLAG(l)  =  1,  one  obtains  these  values  in 
the  file  PROFL.OUT.  Interpolated  terrain  elevation  may  also  be  obtained  in  the  file 
TRNPRF.OUT  if  one  chooses  the  option  IFLAG(2)  =  1 .  I  left  this  option  in  for  future 
graphics  programs.  Currently,  the  screen-graphics  program  GRAPH  creates  contour  plots 
and  by  using  the  option  IFLAG(4)  =  1,  all  points  below  the  terrain  are  given  values  of  500 
dB  losses.  This  way  by  specifying  GL(13)  =  499,  one  obtains  clear  graphical  depiction  of 
the  ground.  Future  codes  will  show  true  losses  and  use  TRNPRF.OUT  to  dravra  in  the 
Earth's  terrain  elevation. 
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Gain  (dB) 


IV  Sample  Input  Files 


The  following  two  tables  contain  samples  of  input  files  used  to  produce  Figures  3  and  4. 
Table  1  gives  the  input  file  for  Example  1,  and  Table  2  gives  the  input  file  for  Example  2.  The 
input  files  can  be  found  on  the  program  CD. 


Table  1:  Example  1  Input  File 

Lloyds  Mirror  TITLE 

4.9965  75.00  75.00  0  0  FREQ  ZS  ZR  POLFLG  ABSRBFLG 
4.000  2.000  4  RMAXDRNDR 

100.00  -100.00  .500  1  100.00  -100.00  HM AX  1  HMAX2DZNDZTPLTBPLT 
30.0  35.0  40.0  45.0  50.0  55.0  60.0  65.0  70.0  75.0  80.0  85.0  90.0 
4  0  110  NPADE  (IFLAG(I),1=1,4) 

3  89.50  ISTRT  THMAX 

.0000  .0000  RD(i)  d(i)  Terrain  Data  Point 

-999  -999 

.0000  .0000  z(i)  N(i) 

-999  -999 

.0000  1.0000  z(i)  EPS(i) 

-999  -999 

.0000  .0000  z(i)  SlGMA(i) 

-999  -999 


Range  (km) 


Figure  3 
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Table  2:  Example  2  Input  File 


Sinusoidal  Terrain 
25.0000  250.00 
50.000  5.000 

1200.00  -500.00 


TITLE 

250.00  1  1  FREQ  ZS  ZR  POLFLG  ABSRBFLG 
10  RMAXDRNDR 

.999  6  1200.00  -200.00  HMAXl  HMAX2  DZ  NDZ  TPLT  BPLT 


60.0  65.0  70.0  75.0  80.0  85.0  90.0  95.0  100.0  105.0  110.0  120.0  499.0 
4  0  0  11  NPADE  (IFLAG(I),1=1,4) 

1  .00  ISTRT  THMAX 

RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d{i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 


.00 

.0000 

.5000 

l.OOOO 

1.5000 
2.0000 

2.5000 
3.0000 

3.5000 
4.0000 

4.5000 
5.0000 

5.5000 
6.0000 

6.5000 
7.0000 

7.5000 
8.0000 

8.5000 
9.0000 

9.5000 
10.0000 

10.5000 
11.0000 

11.5000 
12.0000 

12.5000 
13.0000 

13.5000 
14.0000 

14.5000 
15.0000 

15.5000 
16.0000 

16.5000 
17.0000 

17.5000 
18.0000 

18.5000 
19.0000 

19.5000 
20.0000 

20.5000 
21.0000 

21.5000 
22.0000 

22.5000 
23.0000 

23.5000 
24.0000 


.0000 

61.8034 
117.5571 

161.8034 
190.2113 
200.0000 
190.2113 

161.8034 

117.5570 

61.8034 
.0000 

-61.8034 

117.5571 

161.8034 
190.2113 
200.0000 
190.2113 

161.8034 

117.5570 
-61.8034 

.0000 

61.8034 

117.5571 

161.8034 
190.2113 
200.0000 
190.2113 

161.8034 

117.5570 

61.8034 
-.0001 

-61.8035 

-117.5571 

-161.8034 

-190.2113 

-200.0000 

-190.2113 

-161.8034 

-117.5570 

-61,8033 

.0001 

61.8035 

117.5571 
161.8034 
190.2113 
200.0000 
190.2113 
161.8033 
117.5570 


RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 


24.5000  61.8033 

25.0000  -.0001 

25.5000  -61.8035 
26.0000  -117.5571 

26.5000  -161.8035 
27.0000  -190.2113 

27.5000  -200.0000 
28.0000  -190.2113 

28.5000  -161.8033 
29.0000  -117.5570 

29.5000  -61.8033 

30.0000  .0001 

30.5000  61.8034 
31.0000  117.5571 

31.5000  161.8034 
32.0000  190.2113 

32.5000  200.0000 
33.0000  190.2113 

33.5000  161.8034 
34.0000  117.5570 

34.5000  61.8034 

35.0000  .0000 

35.5000  -61.8034 
36.0000  -117.5571 

36.5000  -161.8034 
37.0000  -190.2113 

37.5000  -200.0000 
38.0000  -190-2113 

38.5000  -161.8034 
39.0000  -117.5570 

39.5000  -61.8034 

40.0000  .0000 

40.5000  61.8034 
41.0000  117.5571 

41.5000  161,8034 
42.0000  190.2113 

42.5000  200.0000 
43.0000  190.2113 

43.5000  161.8034 
44.0000  117.5570 

44.5000  61.8034 

45.0000  -.0001 

45.5000  -61.8035 
46.0000  -117.5571 

46.5000  -161.8034 
47.0000  -190.2113 

47.5000  -200.0000 
48.0000  -190.2113 

48.5000  -161.8034 
49.0000  -117.5570 

49.5000  -61.8033 

50.0000  .0001 

-999  -999 

.0000  376.0000 

200.0000  300.0000 
500.0000  200.0000 
-999  -999 

.0000  5.0000 


RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RX)(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 
RD(i)  d(i)  Terrain  Data  Point 

z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 
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-999  -999 

-500.0000  .0100 
-450.0000  .0010 

-400.0000  .0001 

.0000  .0001 
-999  -999 

5.000000E-01 

61.8034  376.0000 

261.8034  300.0000 

561.8034  200.0000 

-999  -999 

61.8034  5.0000 

-999  -999 

-438.1966  .0100 

-388.1966  .0010 

-338.1966  .0001 

61.8034  .0001 

-999  -999 

1.000000 

117.5571  376.0000 

317.5571  300.0000 

617.5571  200.0000 

-999  -999 

117.5571  5.0000 

-999  -999 

-382.4429  .0100 

-332.4429  .0010 

-282.4429  .0001 

117.5571  .0001 

-999  -999 

1.500000 

161.8034  376.0000 

361.8034  300.0000 

661.8034  200.0000 

-999  -999 

161.8034  5.0000 

-999  -999 

-338.1966  .0100 

-288.1966  .0010 

-238.1966  .0001 

161.8034  .0001 

-999  -999 

2.000000 

190.2113  376.0000 

390.2113  300.0000 

690.2113  200.0000 

-999  -999 

190.2113  5.0000 

-999  -999 

-309.7887  .0100 

-259.7887  .0010 

-209.7887  .0001 

190.2113  .0001 

-999  -999 

2.500000 

200.0000  376.0000 
400.0000  300.0000 
700.0000  200.0000 


z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
zO)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 


Z(i)  EPS(i) 


-999 

200.0000 

-999 


-999 

5.0000 

-999 


-300.0000 

.0100 

z(i)  SIGMA(i) 

-250.0000 

.0010 

z(i)  SlGMA(i) 

-200.0000 

.0001 

z(i)  SlGMA(i) 

200.0000 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

3.000000 

190.2113 

376.0000 

z(i)  N(i) 

390.2113 

300.0000 

z(i)  N(i) 

690.2113 

200.0000 

z(i)  N(i) 

-999 

-999 

190.2113 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-309.7887 

.0100 

z(i)  SIGMA(i) 

-259.7887 

.0010 

z(i)  SlGMA(i) 

-209.7887 

.0001 

z(i)  SIGMA(i) 

190.2113 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

3.500000 

161.8034 

376.0000 

z(i)  N(i) 

361.8034 

300.0000 

z(i)  N(i) 

661.8034 

200.0000 

z(i)  N(i) 

-999 

-999 

161.8034 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-338.1966 

.0100 

z(i)  SIGMA(i) 

-288.1966 

.0010 

z(i)  SlGMA(i) 

-238.1966 

.0001 

z(i)  SlGMA(i) 

161.8034 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

4.000000 

117.5570 

376.0000 

z(i)  N(i) 

317.5570 

300.0000 

z(i)  N(i) 

617.5570 

200.0000 

zO)  N(i) 

-999 

-999 

117.5570 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-382.4430 

.0100 

z(i)  SlGMA(i) 

-332.4430 

.0010 

z(i)  SIGMA(i) 

-282.4430 

.0001 

z(i)  SIGMA(i) 

117.5570 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

4.500000 

61.8034 

376.0000 

z(i)  N(i) 

261.8034 

300.0000 

z(i)  N(i) 

561.8034 

200.0000 

z(i)  N(i) 

-999 

-999 

61.8034 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-438.1966 

.0100 

z(i)  SIGMA(i) 

-388.1966 

.0010 

z(i)  SIGMA(i) 

-338.1966 

.0001 

z(i)  SIGMA(i) 

61.8034 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

5.000000 
.0000  376.0000 


2(i)  N(i) 
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z(i)  N(i) 
z(i)  N(i) 


z(i)  SIGMA(i) 


z(i)  SIGMA(i) 


z(i)  EPS(i) 


-999 

-999 

-438.1967 

.0100 

-388.1967 

.0010 

-338.1967 

.0001 

61.8034 

.0001 

-999 

-999 

15.000000 

-.0001  ; 

376.0000 

200.0000 

300.0000 

499.9999 

200.0000 

-999 

-999 

-.0001 

5.0000 

-999 

-999 

-500.0001 

.0100 

-450.0001 

.0010 

-400.0001 

.0001 

-.0001 

.0001 

-999 

-999 

15.500000 

-61.8035 

376.0000 

138.1965 

300.0000 

438.1965 

200.0000 

-999 

-999 

-61.8035 

5.0000 

-999 

-999 

-561.8035 

.0100 

-511.8035 

.0010 

-461.8035 

.0001 

-61.8035 

.0001 

-999 

-999 

16.000000 

-117.5571 

376.0000 

82.4429 

300.0000 

382.4429 

200.0000 

-999 

-999 

-1  17.5571 

5.0000 

-999 

-999 

-617.5571 

.0100 

-567.5571 

.0010 

-517.5571 

.0001 

-117.5571 

.0001 

-999 

-999 

16.500000 

-161.8034 

376.0000 

38.1966 

300.0000 

338.1966 

200.0000 

-999 

-999 

-161.8034 

5.0000 

-999 

-999 

-661.8035 

.0100 

-611.8035 

.0010 

-561.8035 

.0001 

-161.8034 

.0001 

-999 

-999 

17.000000 

-190.2113 

376.0000 

9.7887 

300.0000 

309.7887 

200.0000 

z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 


2(i)  EPS(i) 


-999  -999 

-190.2113  5.0000 

-999  -999 

-690.2113  .0100  z(i)  SlGMA(i) 

-640.2113  .0010  z(i)  SlGMA(i) 

-590.2113  .0001  z(i)  SlGMA(i) 

-190.2113  .0001  z(i)  SlGMA(i) 

-999  -999 

17.500000 


-200.0000 

376.0000 

z(i)  N(i) 

.0000  : 

500.0000 

z(i)  N(i) 

300.0000 

200.0000 

z(i)  N(i) 

-999 

-999 

-200.0000 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-700.0000 

.0100 

z(i)  SIGMA(i) 

-650.0000 

.0010 

z(i)  SlGMA(i) 

-600.0000 

.0001 

z(i)  SlGMA(i) 

-200.0000 

.0001 

z(i)  SlGMA(i) 

-999 

-999 

18.000000 

-190.2113 

376.0000 

z(i)  N(i) 

9.7887 

300.0000 

z(i)  N(i) 

309.7887 

200.0000 

z(i)  N(i) 

-999 

-999 

-190.2113 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-690.2113 

.0100 

z(i)  SIGMA(i) 

-640.2113 

.0010 

z(i)  SIGMA(i) 

-590.2113 

.0001 

z(i)  SlGMA(i) 

-190.2113 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

18.500000 

-161.8034 

376.0000 

z(i)  N(i) 

38.1966 

300.0000 

z(i)  N(i) 

338.1967 

200.0000 

z(i)  N(i) 

-999 

-999 

-161.8034 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-661.8033 

.0100 

z(i)  SlGMA(i) 

-611.8033 

.0010 

z(i)  SlGMA(i) 

-561.8033 

.0001 

z(i)  SlGMA(i) 

-161.8034 

.0001 

z(i)  SlGMA(i) 

-999 

-999 

19.000000 

-117.5570 

376.0000 

z(i)  N(i) 

82.4430 

300.0000 

z(i)  N(i) 

382.4430 

200.0000 

z(i)  N(i) 

-999 

-999 

-117.5570 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-617.5570 

.0100 

z(i)  SlGMA(i) 

-567.5570 

.0010 

z(i)  SIGMA(i) 

-517.5570 

.0001 

z(i)  SIGMA(i) 

-117.5570 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

19.500000 
-61.8033  376.0000 


z(i)  N(i) 


138.1967  300.0000 

438.1967  200.0000 

-999  -999 

-61.8033  5.0000 

-999  -999 

-561.8033  .0100 

-511.8033  .0010 

-461.8033  .0001 

-61.8033  .0001 

-999  -999 

20.000000 
.0001  376.0000 
200.0001  300.0000 
500.0001  200.0000 
-999  -999 

.0001  5.0000 
-999  -999 

-499.9999  .0100 

-449.9999  .0010 

-399.9999  .0001 

.0001  .0001 
-999  -999 

20.500000 

61.8035  376.0000 

261.8035  300.0000 

561.8035  200.0000 

-999  -999 

61.8035  5.0000 

-999  -999 

-438.1965  .0100 

-388.1965  .0010 

-338.1965  .0001 

61.8035  .0001 

-999  -999 

21.00t)000 

117.5571  376.0000 

317.5571  300.0000 

617.5571  200.0000 

-999  -999 

117.5571  5.0000 

-999  -999 

-382.4429  .0100 

-332.4429  .0010  ' 

-282.4429  .0001 

117.5571  .0001 

-999  -999 

21.500000 

161.8034  376.0000 

361.8034  300.0000 

661.8035  200.0000 

-999  -999 

161.8034  5.0000 

-999  -999 

-338.1966  .0100 

-288.1966  .0010 

-238.1966  .0001 

161.8034  .0001 

-999  -999 


z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SlGMA(i) 
z(i)  SlGMA(i) 
z(i)  SIGMA(i) 


z(i)  N(i) 
z(i)  N(i) 
z(i)  N(i) 

z(i)  EPS(i) 

z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SIGMA(i) 
z(i)  SIGMA(i) 


117.5570 

.0001 

z(i)  SIGMA(i) 

-999 

-999 

24.500000 

61.8033 

376.0000 

z(i)  N(i) 

261.8033 

300.0000 

z(i)  N(i) 

561.8033 

200.0000 

z(i)  N(i) 

-999 

-999 

61.8033 

5.0000 

z(i)  EPS(i) 

-999 

-999 

-438.1967 

.0100 

z(i)  SIGMA(i) 

-388.1967 

.0010 

z(i)  SIGMA(i) 

-338.1967 

.0001 

z(i)  SIGMA(i) 

61.8033 

.0001 

z(i)  SlGMA(i) 

-999  -999 


Range  (km) 


Figure  4 


Good  luck,  and  please  inform  me  of  any  problems  you  may  have.  I  can  be  reached  at 
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Abstract  This  ptq>er  presents  a  solution  of  2D  electromagnetic  wave  propagation  problems 
in  complicated  terrestrial  domains.  Scalar  2D  parabolic  approximations  arc  derived  from 
Maxwell’s  equations  for  both  vertical  and  horizontal  polarization.  The  parabolic  equations 
are  then  solved  using  a  new  technique  involving  Padd  rational  function  approximations  of  the 
macroscopic  operator.  This  method  allows  for  larger-than-normal  range  stepping,  speeding  up 
computational  time  significantly.  The  Padd  approximation  and  the  numerical  implementation 
are  fully  discussed.  The  discontinuity  at  the  earth’s  surface  is  handled  directly  by  using 
classical  continuity  conditions  and  deriving  exact  interface  conditions  for  linking  the  fields 
in  the  atmosphere  to  those  in  the  terrain.  The  interface  conditions  are  then  implemented  using 
the  concept  of  virtual  points.  Preliminary  benchmark  tests  show  the  interface  treatment  to  work 
well.  Finally,  several  example  runs  are  presented  illustrating  results. 


1.  Introduction 

The  solution  of  electromagnetic  (EM)  propagation  problems  in  the  terrestrial  domain  is  a 
complicated  matter.  Three-dimensional  (3D)  variations  in  refraction  and  terrain  make  the  full 
vector  problem  extremely  difficult  to  solve  in  reasonable  time.  If  one  chooses  to  simplify 
the  problem  by  assuming  symmetry  in  one  or  more  of  the  coordinate  directions,  the  vector 
problem  caw  be  uncoupled  into  two  scalar  problems  [1].  However,  the  solution  of  the 
two-dimensional  (2D)  scalar  problem  is  still  very  difficult  for  realistic  environments.  The 
parabolic  approximation  method  is  used  to  reduce  the  solution  of  the  full  two-way  wave 
equation  to  a  solution  of  a  one-way  equation  [2].  Benefits  of  one-way  propagation  ^e  the 
simple  numerical  implementation  of  range  dependencies  in  the  medium,  and  the  avoidance 
of  prohibitive  numerical  aspects  of  solving  elliptic  equations  associated  with  implementing 
two  range-dependent  boundary  conditions.  The  model  discussed  in  this  paper  is  a  so  called 
2.5D  model  using  azimuthally  varied  vertical  planar  fields.  Work  is  also  proceeding  on  a 
3D  model  for  higher  frequencies  using  a  hybrid  combination  of  an  underlying  robust  3D  ray 
trace  and  a  3D  Gaussian  beam  model. 

Two  of  the  most  popular  methods  of  solving  the  parabolic  equations  are  the  implicit 
finite  difference  (IFD)  method  [3]  and  the  split-step  Fourier  transform  method  [4,5].  These 
techniques  are  microscopic  methods  in  the  sense  that  they  are  implementing  approximations 
to  the  differential  equation,  defined  microscopically.  Limitations  of  these  methods  are  that 
the  grid  mesh  over  which  the  solution  is  computed  must  be  small  to  yield  accuracy.  This 
often  means  grid  sizes  on  the  scale  of  at  most  three  wavelengths.  At  high  frequencies  and 
large  propagation  domains  this  could  amount  to  grid  sizes  of  hundreds  of  thousands  by 
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millions  of  points  needed  for  an  accurate  solution.  Methods  using  fourth-order  difference 
schemes  have  also  been  implemented  to  speed  up  computational  time  [6].  A  better  method 
to  solve  the  equations  is  by  symbolically  integrating  the  equation  with  respect  to  range 
thereby  obtaining  the  macroscopic  propagator  [7].  In  theory  if  one  does  an  ‘infinitely’  good 
job  at  approximating  the  macroscopic  propagator  it  is  possible  to  take  infinitely  large  range 
steps.  Limitations  in  adequate  range-dependent  representation  of  the  macroscopic  operator, 
however,  will  limit  actual  range  step  sizes.  In  practice  one  still  can  achieve  great  savings 
in  computational  time  over  IFD  and  even  possibly  the  split-step  Fourier  transform. 

The  macroscopic  operator  is  approximated  using  a  Padd  rational  function  series.  The 
more  terms  used  in  the  approximation,  theoretically  the  larger  the  range-step  one  can  take. 
Another  advantage  of  this  method  is  that  once  the  Pad6  approximation  for  the  propagator 
has  been  obtained  computations  may  be  applied  in  parallel.  That  is,  each  of  the  n  Pad6 
approximations  may  be  applied  to  the  field  at  the  same  time  as  opposed  to  having  to  apply 
a  series  of  products  of  operators  sequentially.  For  cases  when  the  desired  result  is  the  loss 
values  only  near  the  receiver,  and  not  at  many  points  in-between  the  source  and  receiver, 
the  Pad6  technique  applied  to  the  macroscopic  operator  is  ideal.  It  cuts  down  the  number  of 
intermediate  range  locations  at  which  the  field  must  be  computed.  It  is  reiterated,  however, 
that  the  terrain  elevation  and  cover  will  ultimately  determine  the  range  step  size. 

This  report  summarizes  the  theory  involved  in  deriving  parabolic  approximations  to 
scalar  EM  propagation  and  the  associated  boundary  modelling,  including  energy  conservation 
at  vertical  interfaces.  Also  presented  is  the  theory  behind  the  Pad6  rational  function 
approximation  to  the  propagator,  and  the  complete  numerical  implementation  used  in  the 
code  SSP.  We  first  begin  with  the  derivation  of  scalar  wave  equations  for  EM  propagation. 


2.  Derivation  of  scalar  wave  equations 


We  begin  with  Maxwell’s  equations  [1]  in  spherical  coordinates  (r,  0,0),  for  terrestrial 
systems  where  r  is  the  radial  distance  from  the  origin,  6  (measured  positive  down,)  is  the 
angle  between  the  z-axis  and  the  radial  direction  and  0  is  the  azimuthal  angle.  Maintaining 
E  planar  in  the  (r,  6)  plane  leads  to  =  d^£  =  0  for  vertical  polarization  and  H  planar 
in  the  (r,  0)  plane  leads  to  =  d^e  =  0  for  horizontal  polarization.  The  assumption 
of  symmetry  in  the  0  direction  is  a  necessary  requirement  to  reduce  the  original  vector 
problem  to  uncoupled  scalar  problems. 

These  assumptions  used  in  Maxwell’s  equations  readily  lead  to  the  two  equations 


d 

Tdr 


8 

r^sin0  80 


cot0  dn^\ 

Tvle)"* 


=  0 


(Ifl) 


and 
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7^(''^^)+;.2sin0  80 


where  n  is  the  index  of  refraction  defined  as 
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Equation  (la)  determines  the  only  non-zero  component  of  H  in  the  vertical  polarization  case. 
Components  of  E  are  determined  in  terms  of  using  Maxwell’s  equations.  Similarly, 
equation  (li?)  determines  the  non-zero  component  of  E  in  the  horizontal  polarization  case. 
Components  of  H  can  then  be  determined  in  terms  of  using  Maxwell’s  equations.  This 
last  equation  (Ic)  reflects  the  inclusion  of  conductivity  via  current  density  terms  retained  in 
Maxwell’s  equations. 

One  may  transform  equations  (la)  and  (lb)  into  Helmholtz  forms  by  using  the 
substitutions 


H4,= 


nuy 
sin  6 


and 


r  ““ 
Vr  sin  0 


(2a.  b) 


In  doing  so,  equation  (1)  yields 


where 


r  dr  V 


dui 

17 


dh 


302 


-f- 


(*oV- 


4r2  sin^  9 


-|-  Sfii  j  W/  —  0 


Srii  = 


cot0  dn  n  d^n  ’  32n  ’ 

~~7hile  ~  302  ar2 

0 


for  i  =  V 
for  i  =  H. 


(3a) 


(3&) 


We  use  the  subscript  notation  i  =  V(i  =  H)  to  denote  the  vertical  (horizontal)  polarization. 

Since  the  computation  of  rectangular  domains  is  more  desirable  than  spherical  domains, 
we  now  present  an  ‘earth-flattening’  approximation  to  our  problem.  We  will  use  the  smooth 
earth  transformation 


X  =  r^d 


and 


(4a,  b) 


where  is  the  radius  of  the  earth  at  mean  sea  level  (approximately  6370  km).  This 
transformation  places  the  smooth  earth’s  surface  at  z  =  0.  Terrain  will  be  imposed  later 
during  the  numerical  implementation  of  the  solution  algorithm.  The  inverse  transform  is 
given  by 


r  =  Te  exp  ^  ~ 

In  using  the  transformation  defined  in  equation  (4)  differential  operators  translate  as 


and  T - ►  exp 

dr 
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Equation  (3)  thereby  becomes 
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where  m  is  a  modified  index  of  refraction  defined  as 
m  =  ne>.p(i) 
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The  solution  of  equation  (la)  using  parabolic  approximation  techniques  relies  on 
segmenting  the  medium  into  a  series  of  range  independent  sectors.  It  can  be  shown  that  in 
each  of  these  sectors  equation  (la)  becomes 


where 


d^Ui  d^Ui 


+  Kfui  =  0 
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i  =  H 


(8a) 


(89) 


and  the  assumption  of  z/r^  1  is  used.  We  have  also  neglected  the  range-dependent 
cosec^(*)  term  by  assuming  we  are  in  the  far  field  kor  ^  1  and  removed  from  any  poles 
of  the  cosecant  function  (x  <  nr^)  [8].  Equation  (8)  is  the  desired  starting  equation  for 
numerical  implementation.  We  remark  that  modelling  of  media  with  finite  conductivity  is 
achieved  by  replacing  m^  with  the  expression 


+  ifi  ,  (9a 

\CmJ 

where  and  are  modified  light  speed  and  conductivity  defined  as 


Cm  =  cexp(-z/re)  and  (Tm  =  o' exp(2z/re).  (99,  c 


The  next  section  will  discuss  the  boundary  modelling  including  the  effects  anc 
corrections  for  discontinuities  in  range. 


3.  Electromagnetic  boundary  modelling 

Since  we  are  solving  a  parabolic  equation  with  two  z  derivatives  and  one  x  derivative,  wt 
need  two  conditions  in  z  and  one  condition  in  x  at  every  range  step.  Computationally  we  wil 
bound  the  domain  by  two  horizontal  planes  at  the  top  of  the  atmosphere  and  the  botton 
of  the  terrain.  Homogeneous  Dirichlet  conditions  will  be  used  at  these  boundaries  anc 
techniques  will  be  applied  to  reduce  fictitious  reflections.  The  fact  that  the  terrain  introduce; 
a  discontinuity  in  the  problem,  and  its  elevation  is  range  dependent,  complicates  the 
problem  slightly.  Parabolic  approximations  assume  that  range-dependent  environments  art 
partitioned  into  range-independent  sectors.  The  result  is  that  terrain  slopes  are  approximatec 
by  a  series  of  step-like  structures.  This  creates  a  vertical  interface  at  each  range  step  where 
the  terrain  elevation  is  altered  significantly  enough.  Typically,  the  condition  at  vertica 
interfaces  is  continuity  of  field,  however  this  is  not  the  best  condition  one  might  impose 
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We  will  eventually  impose  a  conservation  of  energy  condition  at  discontinuities  in  range  to 
correct  for  some  range-dependent  errors. 

We  first  discuss  the  modelling  of  the  horizontal  interface  between  the  atmosphere  and  the 
earth’s  surface.  Assuming  finite  conductivity  and  no  surface  charge  density  the  conditions 
on  B,  H,  D  and  E,  at  an  interface  between  two  different  media  (Sa,  and  (Sb.  Mb.  c^b) 

are  given  in  [1].  In  the  case  of  a  smooth  interface,  the  normal  vector  e„  =  e,.  Assuming 
also  that  magnetic  permeability  in  the  terrain  is  that  of  the  atmosphere  (both  equal  to  the 
free  space  value,)  these  conditions  imply  the  following  six  relations  hold: 

(10fti.,c) 

Hea  =  Heb  and  Hru  =  Hrb-  (10^.  ^J) 

For  efficient  numerical  computation,  the  approach  here  treats  the  variable  terrain  as 
a  series  of  discrete  vertical  jumps,  i.e.  a  staircase  approximation,  using  smooth  earth 
formulae  for  the  horizontal  boundary  conditions  in  each  sector.  However,  the  general 
boundary  conditions  for  the  variable  terrain  have  been  determined  for  vertical  and  horizontal 
polarization.  The  implementation  of  the  general  conditions  is  being  examined. 

For  the  case  of  vertical  polarization,  clearly  one  condition  on  is  given  by 
equation  (lOrf).  The  second  condition  may  be  obtained  by  manipulating  Maxwell’s  equations 
and  equation  (lOi?)  yielding 


£a  dr  Bb  or 

at  the  terrain  surface.  A  third  condition  is  also  given  by  dH^JdO  =  dH^b/d6  although 
this  will  not  be  used.  For  the  case  of  horizontal  polarization,  equation  (10a)  will  be 
one  condition,  and  the  second  is  obtained  by  manipulating  Maxwell’s  equations  and 
equation  (lOe)  giving 

■^{rE^a)  =  ■|-(r£^b) 
or  or 

at  the  terrain  surface.  The  reason  that  the  horizontal  polarization  case  ‘does  not  see’  a 
discontinuity  is  because  there  is,  in  fact,  no  discontinuity  in  permeability  at  the  surface. 

We  next  use  the  flattening  transformation  defined  in  equations  (4),  and  equations  (11) 
and  (12)  to  give 
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We  now  switch  independent  variables  using  equation  (2).  In  the  case  of  vertical 
polarization  one  derives  from  equation  (lOJ) 


«aMVa  =  ”b“Vb 


(14a) 
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while  equation  (13a)  gives 
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Switching  to  modified  refiraction  m,  in  using  equation  (lb)  these  equations  become 

mallya  =  //IbMvb  (15a) 


and 
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For  the  case  of  horizontal  polarization,  use  of  equation  (2i>)  and  equations  (10a)  and 
(12)  gives 


MHa  =  MHb  and 


9MHa  _  9MHb 
9z  9z 


(16a.  6) 


For  efficient  numerical  computation,  the  approach  taken  here  treats  the  variable  terrain  in 
terms  of  staircase  approximations  using  smooth  earth  formulae  for  the  horizontal  boundary 
conditions.  The  smooth  earth  boundary  conditions  are  simplifications  of  the  general 
boundary  conditions  for  variable  terrain.  These  conditions  can  be  determined  without 
approximations  to  the  terrain,  for  the  cases  of  vertical  and  horizontal  polarization,  in  terms 
of  directions  normal  and  parallel  to  the  terrain.  Using  a  generalization  of  the  earth  flattening 
transformation  that  is  applicable  to  the  variable  terrain,  the  general  conditions  can  also  be 
given  in  earth  flattened  coordinates  {x,  z).  For  brevity  these  formulae  are  not  included. 

We  now  discuss  the  discontinuity  at  vertical  interfaces,  whether  real  (due  to  range 
dependent  changes  in  refraction)  or  numerical  (due  to  stair  case  approximations  to  terrain 
slopes.)  Any  parabolic  equation  requires  two  conditions  in  height  and  one  condition  in 
range  to  determine  a  unique  solution.  We  have  just  derived  the  two  height  conditions  at 
horizontal  interfaces  for  each  polarization  case.  When  considering  the  condition  at  vertical 
interfaces,  once  an  initial  field  has  been  specified,  at  cacti  step  in  range  the  typical  condition 
is  continuity  of  field.  That  is  to  say  that  when  a  .'ct  of  loss  values  has  been  computed  at 
a  particular  range  step,  these  values  are  used  as  initial  data  for  taking  the  next  range  step. 
This  condition  can  be  replaced  by  other  possible  more  desirable  conditions.  We  use  results 
from  scalar  acoustic  problems  to  derive  a  conservation  of  energy  condition  for  the  vertical 
interface  [9-12]. 

The  basic  concept  is  to  introduce  a  new  dependent  variable  obtained  by  scaling  the 
old  variable.  The  scale  could  be  height-dependent.  In  order  to  properly  implement  an 
energy-conserving  condition  at  a  vertical  interface  we  must  derive  an  expression  for  the 
energy  flux  in  the  range  direction.  We  begin  with  the  complex  Poynting  vector  S,  which 
is  generally  taken  to  give  the  flow  of  energy  in  a  propagating  electromagnetic  field.  The 
vector  is  defined  as  [1] 


S  =  \ExH*  (17) 

where  the  asterisk  denotes  complex  conjugate.  The  average  intensity  of  energy  flow  is  taken 
as  the  real  part  of  the  complex  Poynting  vector.  For  the  case  of  vertical  polarization 

S={{{EeH;)er-{ErH;)ee). 


(18) 


(19) 


The  energy  flux  in  the  B  direction,  denoted  Sg,  is 
Using  Maxwell’s  equations  and  equation  (2)  in  equation  (19)  gives 


5fl  = 
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Im 


'duv  *■ 


2coEor'^sin6  V 


(20) 


In  using  equation  (4)  to  switch  to  earth  flattened  coordinates  the  energy  flux  is  given  by 

The  condition  for  conservation  of  energy  along  a  vertical  interface  thereby  becomes 


Im 


.  dx 


(22) 


where  we  use  the  subscript  ‘in’  here  to  denote  the  incident  wave  while  ‘tr’  denotes  the 
transmitted  wave. 

This  type  of  nonlinear  boundary  condition  is  quite  difficult  to  implement  in  practice, 
however  an  equivalent  linear  condition  can  be  derived.  By  factoring  equation  (8a),  and 
retaining  the  outgoing  solution,  one  can  derive 


9«v 

dx 


=  ikoQy/U\ 


(23a) 


where 


koQw  =  Y  ^  +  -^v- 
Using  equation  (23)  in  equation  (22)  gives 

Im(il2Vtfl«Vul^)  =  I>^(il2vi„lwVi„l^) 


(23b) 


(24a) 


which  will  be  satisfied  if 


1/2, 


(24b) 


One  could  in  theory  apply  equation  (24b)  as  the  propagator  for  steps  across  vertical 
discontinuities.  The  simpler  method  of  incorporating  conservation  conditions  via  a  new 
dependent  variable  is  derived  by  assuming  negligible  propagation  angles.  With  this 
assumption,  and  those  to  follow,  equations  (8fe)  and  (23b)  yield  the  approximation  Q  « 
(Ky/kiy^^  «  Vm  so  that  the  condition  becomes 

=  V^MVi. 


where  m  is  a  modified  index  of  refraction  defined  by  equation  (9).  In  deriving  equation  (25), 
one  must  assume  that  the  terms  involving  z  derivatives  of  m  ‘  are  negligible.  We  obtain  here 
merely  a  first-order  correction  for  conserving  energy.  If  one  desired,  one  could  implement 


the  full  correction,  given  in  equation  {2Ab)  at  each  range  step  where  there  is  a  vertical 
discontinuity  of  some  type.  This  compounds  the  numerical  process,  however,  and  slows 
down  the  algorithm  significantly.  It  has  been  demonstrated  that  the  first-order  correction  is 
sufficient  for  many  types  of  environments  when  considering  acoustic  propagation  [9].  Tests 
underway  also  tend  to  support  this  conclusion  for  electromagnetic  problems  as  well. 

Equation  (25)  suggests  the  transformation  of  dependent  variable  to  be 

xuy  =  -/muy  (26) 

and  continuity  of  field  in  w  will  imply  a  first-order  correction  for  conservation  of  energy  in 
u.  The  case  of  horizontal  polarization  is  analysed  in  a  similar  fashion  leading  to  analogous 
results  and  a  transformatica  identical  to  equation  (25).  For  brevity,  these  results  are  omitted. 

In  using  this  transformation  for  conserving  energy  the  elliptic  equation  that  is  solved  is 
given  as 


where  this  equation  has  been  derived  by  substituting  equation  (26)  into  equation  (8a).  In 
implementing  the  energy-conserving  transformation  one  must  now  transform  the  interface 
conditions  given  in  equations  (15)  and  (16).  For  the  case  of  vertical  polarization  conditions 
on  wy  become 

•v/^iuva  =  -v/^uivb  (28a) 


and 
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The  conditions  on  ioh  fo;  the  case  of  horizontal  polarization  become 
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and 
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The  preceding  derivation  has  assumed  that  the  complex  part  of  m  is  small  compared  to 
the  real  part.  This  is  true  over  most  of  the  frequency  regime  and  terrain  cover  of  interest. 


4.  Numerical  implementation 

Our  goal  is  to  solve  the  differential  equation,  equation  (27), 


subject  to  homogeneous  Dirichlet  conditions  at  the  boundaries,  and  the  interface  conditi^s 
in  equation  (28)  for  vertical  polarization  or  equation  (29)  for  horizontal  polarization.  The 
modified  wavenumber.  is  defined  in  equation  (8t).  For  convenience  in  u^oming 
notation,  we  drop  the  i  subscript  in  favour  of  the  reader  understanding  that  all  results 
apply  to  both  the  vertical  and  horizontal  polarization  cases.  This  equation  may  be  formally 

factored  to  give 


^  =  i*„yiT2^ 

dx 


(30a) 


where 


(306) 


as  the  differential  equation  governing  the  outwardly  propagating  wave.  The  differences 
between  this  formulation  and  that  in  equation  (23)  is  numerically  motivated  and  will  be 
discussed  shortly.  Removing  the  exp(ikox)  from  the  solution  u,  equaUon  (30a)  becomes 


—  =  iito(-l  +  \/l  +  l2)j"- 
aj:  V  '' 

One  may  formally  integrate  equation  (31)  to  give 

w(x  +  Ax,  z)  =  exp  ^i/:oAx  ^-1  +  \/r+  w(x,  z). 
Following  the  method  of  Collins  we  apply  a  Pad6  approximation  [7]. 

_  np  ^ 

exp  (ikoAx  (-1  +  ^/lTe))  =  1  +  E  i  +  l[„pQ 


(31) 


(32) 


(33) 


where  the  coefficients  ai,„p  and  bi,np  are  determined  numerically  using  the  approach  in  [11]. 
The  method  used  in  [11]  converges  faster  using  the  formulation  in  equation  (30)  than  that 
in  equation  (23).  The  number  np  is  the  Pad6  number,  or  the  number  of  terms  used  in  the 

series  approximation.  .  r.  i  • 

Substituting  equation  (33)  into  equation  (32)  one  obtains  the  split-step  Pade  solution, 


w(x 


-f  Ax,z)  =  w(x,z)  + 


fl/.np  z) 

^  1  +  bi^np  Q 


(34) 


The  terms  in  the  sum  may  be  computed  in  parallel,  which  is  what  makes  this  technique  so 
appealing.  That  is  we  compute 


fi(x,z) 


ai,npQMx,z) 

1  +  bi^np  Q 


(35a) 


in  parallel  and  then  calculate 

np 

w(x  -t-  Ax,  z)  =  w(x,z)  +  y~l  ^i(x,  z). 

/=! 


(356) 
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Let  us  discretize  the  problem  as  follows.  First  we  will  use  a  simple  linear  transformation 
to  invert  the  problem  so  that  z  is  measured  positive  down  from  the  top  of  the  atmosphere. 
For  ease  of  notation  we  will  still  use  the  variable  z  rather  than  defining  a  new  variable,  say 
z'.  Then  we  define  a  grid  with  mesh  sizes  dx,  and  dz.  Let 

w"  —  windx,  j dz)  (36a) 

and 

^i.j  =  dz)  (36^) 

■yVnen  discretizing  the  differential  operator  Q,  equation  (35a)  becomes  a  tridiagonal 
linea:  system  of  which  the  7th  equation  is 

Rlijylflj.,  +  Rlijflj  +  R\j^lj+x  =  +  Shjw’l  +  53,,yu;;+,  (37) 

where  /?1,  Rl,  R3,  51,  52,  and  53,  are  dependent  upon  dz  and  medium  properties  through 
the  function  K  and  m.  Once  this  system  is  solved  for  xfrJ^j  one  uses  equation  (35b)  to 

compute  the  solution  w"'^^  as 


np 

/=i 

The  numerical  domain  is  terminated  with  homogeneous  Dirichlet  conditions  on  the 
field  w.  To  avoid  spurious  reflections  from  the  top  of  the  atmosphere,  an  absorbing  layer 
is  introduced  with  complex  wavenumber.  Similarly,  at  the  bottom  of  the  earth  layer,  one 
increases  conductivity  so  as  to  eliminate  reflections.  For  most  typical  ground  cover  with 
non-zero  conductivity  the  Earth  acts  as  an  absorbing  layer  naturally.  Increasing  conductivity 
near  the  bottom  of  the  domain  insures  no  reflections.  As  for  implementing  the  interface 
conditions  as  given  in  equations  (28)  and  (29)  one  uses  the  idea  of  virtual  points  [11). 
Assume  the  atmosphere-terrain  interface  occurs  between  the  jth  and  (y-f  l)th  notch  points. 
We  n*ace  two  virtual  points  a!  and  b'  in  between  the  two  actual  points,  such  that  the  point 
a'(b')  represents  the  continuation  of  the  atmospheric  (terrain)  solution  one  notch  point. 

This  technique  is  described  fully  in  [10].  By  requiring  each  to  satisfy  the 
linear  interface  conditions  at  the  earth’s  surface,  one  automatically  satisfies  the  interface 
requirement  on  the  entire  solution  lu".  The  discretized  equations  at  the  nodes  on  each  side 
of  the  interface  are 


Rhjirl'j.,  -f  Rlijfl'j  +  R3ijtl^  =  +  53, 
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We  approximate 
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at  the  interface,  and  substituting  these  expressions  in  equation  (28)  for  vertical  polarization 
or  equation  (29)  for  horizontal  polarization  allows  one  to  solve  for  the  virtual  point  solutions 
in  terms  of  actual  notch  point  solutions.  For  the  case  of  vertical  polarization,  using 
equation  (40)  in  equation  (28)  gives 
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Solving  equations  (41a)  and  (41&)  for  the  virtual  solutions  gives 
w",  =  aijWj  +ai2Wj^i 


and 


=a2iwj +a22w"_^j 


where 

ail  =  (AifAs  —  Ae)  +  A2(A4  —  A3))/(Ai(A6  -  As)  +  A2(A3  +  A4)) 
ai2  =  (2A2A6)/(Ai(A6  —  As)  +  A2(A3  +  A4)) 

“21  =  (2Ai  A4)/(Ai(A6  —  As)  +  A2(A3  +  A4)) 
and 

“22  =  (Ai(As  +  Ab)  —  A2(A3  +  A4))/(Ai(A6  —  As)  +  A2(A3  +  A4)). 


(41c,  d) 
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One  may  also  repeat  this  procedure  for  each  individual  since  they  each  satisfy 

the  same  interface  conditions  one  uses  equation  (42)  exactly  for  virtual  point  solutions 
and  The  result  is  that  at  the  interface  equation  (39)  becomes 

+  (R'^i.j  +  otnR\j)flj  +ocnR'iijf”j+i 

=  ShjwU  +  iS2ij  +  ccuS3ij)w]  +  anS\jw'l+^ ,  (43a) 


and 

a2\  Rli,j+\xlf"j  +  (R2ij+\  +  oi22Rh,j+\W"j+i  +  R^ij+xi^ij+i 

=  a2iSlij+\u}j  +  (S2ij+\  +“2251i./+i)ia"+i  +  •53/j+iw"+2* 


(43h) 
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A  similar  analysis  in  the  horizontal  polarization  case  yields  very  similar  results  in 
terms  of  equation  (43),  with  differences  in  the  definitions  of  the  a, 7  coefficients.  Using 
equation  (40)  in  equation  (29)  gives  the  equations 


A\(Wj  +  u)"/)  —  A2(iw”+j  4- 

(44a) 

and 

+  0  +  Ai(iyJ  -  tup  =  ^(u>p,  +  <)  +  A2(wp, 

-  (446) 

where  now 
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Solving  equations  (44a)  and  (44fo)  for  wl  and  uij,  yields 

=  «1I 

(45a) 

and 

lOb-  =  Ci2\  Wj  +  a22w]+} 

(456) 

where 

ail  =  ^22  =  (A1A4  —  A2A3)/(2Ai  A2  +  A2A3  —  Ai  A4) 

(45c) 

ai2  =  (2A|)/(2AiA2  +  A2A3  -  Ai  A4) 

(45d) 

and 

a2i  =  (2A;)/(2Ai  A2  +  A2A3  -  Ai^). 

(45c) 

The  final  implementation  is  exactly  that  in  equation  (43)  with  or/y  being  replaced  by  a, 7. 

Therefore  the  numerical  implementation  of  the  terrain  interface  is  a  simple  modification 
of  the  algorithm  at  the  ;th  and  {j  +  l)th  nodes.  The  code  SSP  is  currenUy  undergoing 
testing.  The  next  section  discusses  preliminary  testing  and  evaluation  of  SSP  including  the 
atmosphere-terrain  interface  implementation. 


5.  Numerical  examples 

The  numerical  code  SSP  is  currently  being  tested.  The  results  presented  in  this  paper  do  not 
utilize  the  advantage  of  parallel  processing.  Preliminary  solutions  were  calculated  on  PCs 
simply  to  test  the  code  and  demonstrate  certain  aspects  of  this  method.  The  final  program 
will  be  run  on  a  parallel  machine.  References  [8-10]  suggest  run  time  speeds  100  times 
faster  than  conventional  methods. 
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Since  the  implementation  of  interface  matching  conditions  is  an  integral  part  of  the 
program,  examination  of  the  numerical  methods  used  there  is  critical.  While  some  acoustic 
propagation  models  do  incorporate  horizontal  interfaces  many  more  do  not.  It  seems  that 
the  discontinuity  in  the  acoustic  case  is  slight  as  compared  to  the  EM  case.  So  slight  as  to 
permit  one  to  smooth  out  discontinuities  or  simply  ignore  them  altogether.  This  is  certainly 
not  the  case  for  electromagnetic  propagation.  A  series  of  tests,  too  comprehensive  to  present 
here,  have  examined  this  aspect  of  the  problem.  By  altering  acoustic  parameters  so  as  to 
have  a  discontinuity  of  order  similar  to  the  EM  case,  the  results  show  that  one  must  also 
incorporate  the  acoustic  matching  conditions  to  have  accurate  solutions. 

The  simplest  analytical  model  to  test  the  method  is  the  Lloyds  mirror  probleni,  in  which 
th6  medium  is  taken  as  homogeneous  with  homogeneous  Dirichlet  conditions  at  the  top 
and  bottom  of  the  domain.  Tlie  source  frequency  is  taken  as  4.99  MHz,  with  dr  =  2  m, 
and  dz  =  0.5  m.  The  source  and  receiver  height  are  both  taken  to  be  75  m,  and  the  entire 
thickness  of  the  domain  is  200  m.  Figure  1(a)  shows  a  comparison  of  the  propagation 
losses  as  computed  by  the  program  SSP  (full  curve)  and  the  program  efepe  (dotted).  The 
code  EFEPE  is  the  original  acoustic  code  and  was  benchmarked  against  a  normal  mode 
program  showing  excellent  agreement  [7].  As  one  can  see  there  is  virtually  no  difference 
in  solutions.  We  were  satisfied  that  the  EM  adaptation  was  coded  correctly.  We  next 
considered  the  interface  modelling.  We  are  currently  looking  into  benchmark  models  for 
EM  propagation  so  as  to  fully  test  the  numerical  implementation  of  the  interface  matching 
conditions.  However,  there  are  existing  acoustic  benchmark  models  available  immediately. 

We  are  able  to  implement  acoustic  boundary  conditions,  similar  to  our  electromagnetic 
conditions,  in  the  code  and  compare  them  to  results  that  were  benchmarked  against  normal 
mode  solutions.  The  model  we  choose  is  called  Case  3b  from  the  NORDA  Parabolic  Equation 
Workshop  [13].  It  is  basically  one  homogeneous  medium  overlying  a  different  homogeneous 
medium  with  absorption.  The  density  discontinuity  at  the  ocean-sea  floor  interface  is  very 
similar  to  the  light  speed  discontinuity  in  the  EM  problem  at  the  atmosphere-terrain  interface. 
While  the  interface  conditions  for  the  acoustic  model  are  not  as  complicated  as  the  EM 
problem,  the  numerical  implementation  using  virtual  points  is  identical. 

Having  implemented  the  interface  conditions  in  the  acoustic  code,  figure  1(b)  compares 
the  results  from  SSP  (full)  and  the  oiiginal  acoustic  program  EFEPE  (dotted,)  which  was 
benchmarked  against  a  normal  mode  solution.  The  output  from  EFEPE  was  in  excellent 
agreement  with  the  normal  mode  solution  except  at  null  locations,  specifically  the  one  near 
7  km.  While  there  are  subtle  differences  in  losses,  they  are  at  most  1  dB,  and  the  curves 
generally  agree  quite  well.  Work  is  currently  underway  to  test  the  EM  interface  conditions. 

The  remaining  figures  demonstrate  the  method’s  ability  to  greatly  increase  range  step 
size.  We  will  use  the  Lloyds  mirror  example  to  illustrate.  Using  the  EFEPE  result  in 
figure  1(a)  as  a  benchmark  we  have  calculated  the  solution  using  the  code  SSP  with  np  =  4 
and  dr  =  50  m.  Figure  2(a)  shows  the  results  with  fair  agreement  that  decays  as  range 
increases.  When  increasing  np  to  8,  as  in  figure  2(b),  excellent  accuracy  is  obtained. 

For  parallel  computations  this  could  mean  very  little  extra  run  time  achieving  much  better 
accuracy.  When  the  range  step  is  increased  to  100  m,  as  in  figure  3(a),  the  accuracy  is  still 
maintained  for  np  =  8.  When  the  range  step  is  increased  to  200  m,  the  accuracy  begins 
to  degrade.  Finally,  when  increasing  np  to  10,  shown  in  figure  4,  accuracy,  while  still  not 
perfect,  is  increased  greatly. 

In  theory,  one  can  take  very  large  range  steps  when  using  this  method.  One  need  only 
take  np  large  enough  so  that  the  Pad6  approximation  gets  arbitrarily  accurate.  However,  the 
limiting  factor  will  be  terrain  and  atmospheric  conditions.  Taking  range  steps  too  large  could 
result  in  ‘stepping  over  mountains’.  As  is  typical  with  parabolic  approximation  methods, 
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Figure  1.  (a)  The  Lloyds  mirror  test  problem.  Source 
and  receiver  height  are  75  m,  and  source  frequency  is 
4.99  MHz.  Intensity  loss  in  dB  is  plotted  against  range 
in  km.  ssp  (full  curve)  and  efepe  (dotted  curve)  show 
excellent  agreement,  (b)  Acoustic  benchmark  Case  3b 
in  [13]:  ssp  (full  curve)  and  efepe  (dotted  curve)  show 
excellent  agreement. 


Figure  2.  Lloyds  mirror  test  problem:  comparison 
of  SSP  with  benchmark  solution  for  (a)  np  =  4  and 
dr  =  50  m,  and  (b)  np  =  8,  and  dr  =  50  m. 


the  range-dependent  problem  is  sectored  into  range-independent  slabs.  Discretization  of 
a  continuous  medium  results  in  staircase  effects.  These  effects  are  minimized  by  taking 
smaller  range  steps.  Therefore  there  is  still  a  lot  of  work  to  be  done  in  the  delicate  matter 
of  trading  off  time  (dr  and  np  as  large  as  possible)  and  accuracy  (dr  small  enough  to 
capture  the  true  physical  properties  of  the  medium).  Sensitivity  testing  and  benchmarking 
are  crucial  aspects  of  this  problem  and  current  efforts  are  being  placed  on  these  areas. 


6«  Summary 

Scalar  Helmholtz  equations  have  been  derived  directly  from  Maxwell’s  equations  for 
the  cases  of  vertical  and  horizontal  polarization.  The  primary  assumption  necessary  for 
such  a  reduction  is  that  the  medium  is  approximately  symmetric  in  one  spatial  direction. 
Factorization  of  these  equations  yield  parabolic  equations  which  are  then  symbolically 
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Figure  3.  Lloyds  mirror  test  problem:  comparison 
of  SSP  with  benchmark  solution  for  (a)  n/7  =  8  and 
dr  =  100  m,  and  (b)  np  =  8,  and  dr  =  200  m. 


Figure  4.  Lloyds  mirror  test  problem:  comparison 
of  SSP  with  benchmark  solution  for  np  =  10  and 
dr  =  200  m. 


solved  in  the  range  direction.  What  one  obtains  is  a  symbolic  expression  for  the  range¬ 
stepping  macroscopic  operator.  Rather  than  discretize  the  microscopic  operator,  a  Pad6 
series  approximation  is  used  for  the  macroscopic  operator.  In  theory  this  allows  very  large 
range  steps.  Step  size  is  still  ultimately  determined  by  medium  characteristics.  However 
there  are  numerical  advantages  over  typical  finite  difference  methods,  the  main  being  the 
suitability  of  the  method  to  be  parallelized  for  multi-processor  computers. 

Interface  conditions  have  been  fully  developed  for  linking  the  atmosphere  to  the 
ground.  Methods  for  conserving  energy  at  vertical  interfaces  has  also  been  discussed, 
with  the  result  that  to  first  order,  a  simple  transformation  of  dependent  variable  allows 
for  implementation.  This  will  provide  corrections  to  sloping  terrain  errors  as  well  as 
range-dependent  refractive  effects.  The  numerical  implementation  of  the  split-step  Pad6 
solution  and  interface  conditions  has  also  been  presented.  Several  benchmark  calculations 
and  interface  modelling  comparisons  were  also  presented.  A  full  description  of  the  code 
SSP  is  contained  in  technical  reports  available  from  the  authors  upon  request.  The  user’s 
manual  includes  program  flow  charts,  input  descriptions  and  output  options.  Also  described 
in  this  report  is  a  post  processing  graphics  program  called  GRAPH.  This  program  produces 
contour  graphs  for  visual  display  only. 
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Appendix  B 

Absorbing  Layer  Validation 


Absorbing  Layers 


A  major  complication  in  the  solution  of  propagation  problems  in  physically  infinite 
media  is  the  attenuation  of  fictitious  reflections  from  computational  boundaries  necessary 
to  terminate  the  calculation.  For  fuU  elliptic  or  hyperbolic  problems  the  boundaries  are  not 
only  the  top  or  bottom  of  the  domain  but  the  right-hand  edge.  For  a  parabolic  problem  one 
need  only  worry  about  the  top  and  bottom  edges  There  are  two  approaches  to  solving  such 
problems.  The  first  technique  is  to  use  an  impedance-type  boundary  condition  at  the 
computational  edge.  The  second  technique  places  an  absorbing  layer  at  the  region  in 
question.  This  layer  has  the  property  that  the  index  of  refraction  has  an  imaginary 
component  that  varies  with  height.  The  absorbing  potential  (the  imaginary  piece)  must  vary 
in  such  a  way  as  to  eliminate  reflections.  The  function  has  a  maximum  at  the  end  of  the 
absorbing  layer  (in  our  case  near  the  top  of  the  domain)  and  decreases  exponentially  as  the 
distance  from  the  boundary  increases.  I  have  obtained  the  foUowing  information  from  Dr. 
Steven  Wales,  at  N.R.L.,  and  I  (he)  has  no  explanation  as  to  why  it  works  but  after  years  of 
experience  this  is  what  he  and  NRL  came  up  with.  The  imaginary  part  of  the  index  of 
refraction  is  taken  as 


Vo  exp 


(i6xy  j 


where  is  taken  to  be  .10  (in  mks  system).,  and  X  is  source  wavelength.  In  terms  of  the 
reference  wave  number  ,  the  fiinction  is 


Vq exp(-9. 8946468 X  10‘*^o(^~^B  )^) 


Dr.  Wales  had  suggested  that  for  most  frequencies  of  interest  the  width  of  the 
absorbing  layer  should  be  about  1/3  the  thickness  of  the  actual  medium,  however,  he  also 
mentioned  that  at  the  lower  frequencies  the  width  had  to  be  greater.  I  began  testing  the 
layer  by  picking  various  width  factors.  After  several  test  runs  I  have  determined 
appropriate  width  factors  for  frequencies  of  1,  3,  4,  5,  10,  20,  and  25  MHz.  It  seems  that 
for  frequencies  above  25  MHz,  the  choice  of  a  width  factor  =  1/3  is  in  fact  good.  However, 
for  the  lower  frequencies  several  values  had  to  be  selected.  The  results  are  summarized  in 
the  following  table. 


Frequency 

(MHz) 

Wave  number 

(^o) 

Width  Factor 

1 

.020958 

15 

3 

.062874 

6 

4 

.083832 

4 

5 

.010479 

2 

10 

.20958 

1 

20 

.41916 

.56 

25 

.52396 

1/3 

By  width  factor,  I  mean  the  amount  of  the  requested  physical  domain  that  is  extended.  If 
HMAXl  is  750,  and  the  width  factor  is  1/3  then  the  absorbing  layer  is  computed  as  250  ir 
wide.  If  HMAXl  is  1000m,  and  the  width  factor  is  13,  then  the  absorbing  layer  is  taken  tc 
be  13,000m.  The  following  pages  show  comparisons  of  the  code  SSP  with  the  analytica 
solution  for  an  infinite  homogeneous  half  space  Dirichlet  problem  .  The  solution  to  the 
spherically  symmetric  wave  equation  is 


u  =  -e\p{ikop) 

P 

where  p  is  the  distance  fi’om  the  source  to  the  receiver. 

Ignoring  curvature  effects,  let  a  cylindrical  coordinate  system  be  placed  so  that  . 
source  is  positioned  on  the  z-axis  at  (r,Zj).  If  the  earth  (at  z  =  0,)  is  such  that  t/  =  0  there 
one  may  use  the  method  of  images  to  detemune  that 


where 


u  = 


— exp(/  ^0  A  )  -  —  exp  (/  ko  P2 
Pi  P2 


). 


Pi  =ylr^  +(z-z,y. 


and 


A  =  V'"'  +  • 


Figure  1-7  show  a  comparison  of  SSP  and  the  analytical  solution  for  frequencies  of 

l,  3,  4,  5,  10.  20,  and  25  MHz.  In  each  case  there  are  three  figures  ( i.e.  Fig  (la),  (lb)  and 
(Ic)  ),  Figure  #(a),  shows  the  entire  interval  0  <  r  <  50  km,  while  Fig.  #b,  and  #c  show 
blow  ups  of  the  range  intervals  0  <  r  <  2  km  and  48  <  r  <  50  km  respectively.  Clearly 
there  is  excellent  agreement.  The  source  and  receivers  are  positioned  250  m  above  the 
surface.  I  ch  -ose  the  medium  to  be  2000  m  high.  There  are  some  phase  errors  early  on,  but 
I  cannot  be  sure  if  the  errors  are  due  to  the  absorbing  layer  or  the  inherent  problems  that 
PE's  have  in  the  near-field.  In  either  case  the  phase  errors  disappear  soon  afterwards.  For  a 
source  frequency  of  1  MHz,  as  seen  in  Fig  (Ic),  there  are  differences  of  at  most  O(.Ol)  dB. 
This  is  also  true  for  /=  3  MHz  (Fig.  (2c))  and  /  =  4  MHz  (Fig.  (3c)).  For  source 
frequencies  of  5,  10  and  20  MHz,  the  differences  are  even  less  than  the  preceding  cases. 
These  frequencies  correspond  to  acoustic  frequencies  of  approximately  25,  50,  and  100  Hz. 
This  was  the  frequency  range  of  interest,  I  believe,  to  NRL  and  Steven  Wales.  For  higher 
frequencies,  the  agreement  begins  to  decay  slightly  as  seen  in  Fig.  (7c). 

To  check  the  programs  interpolation  of  layer  widths,  I  ran  test  cases  for  frequencies 
of  2,  3.5,  4.5,  7.5,  15,  and  22.5  MHz.  These  results  are  shown  in  Fig.  8  -  14,  with  similar 
range  intervals  as  in  the  first  7  graphs.  The  results  show  exceUent  agreement  at  the  longer 
ranges  with  dB  differences  typically  of  the  order  .01.  In  some  cases  there  are  differences  of 
up  to  .15  dB,  but  that  is  still  quite  excellent. 

Figure  14  shows  results  for  a  source  frequency  of  50  MHz.  This  is  the  case  where  I 
saw  the  greatest  difference.  Clearly,  there  should  be  some  modification  of  the  ampUtude 
and  half-width  of  the  Gaussian  absorbing  potential  as  frequency  is  increased.  However, .  1 5 
dB  is  raiher  a  small  difference  at  25  km  which  is  more  than  4000  wavelengths.  For  fiature 
work,  I  suggest  that  someone  does  look  into  modification  of  the  absorbing  potential  for 
higher  frequencies. 

I  also  ran  a  few  cases  with  different  source  and  receiver  locations  for  long  ranges. 
Figure  15  show  results  for  a  source  frequency  of  6.2  MHz.  The  source  is  located  at  182  m 
above  the  surface  and  the  receiver  is  262  m  high.  The  maximum  height  (HMAXl)  is  3000 

m.  As  seen  in  Fig  (15c),  the  solution  begins  to  decay  near  r  =  100  km.  It  is  clear  that  while 
the  differences  are  less  than  1  dB,  the  oscillations  are  increasing.  The  differences  are 
resolved  by  raising  HMAXl  to  5000  m,  as  in  Fig  (15d),  and  errors  are  now  again  between 
.01  and  .  1  dB.  Figure  16  is  a  situation  where  the  source  frequency  is  72. 1  MHz  and  source 
height  is  151  m.  The  receiver  is  located  at  82  m  above  the  surface  and  HMAXl  =  2000  m. 
While  the  maximum  range  is  50  KM  not  100  km  as  in  the  previous  example,  the  number  of 


wavelengths  is  approximately  12,500.  This  is  a  very  large  problem.  As  seen  in  Fig  (16c), 
there  is  again  excellent  agreement  between  the  analytic  and  computed  solutions. 

I  ran  these  odd  number  cases  just  to  make  sure  I  wasn't  getting  any  fimny  results  b 
placing  the  source  and  receiver  at  250  m  which  was  at  a  depth  node  (harmonically  spearir 
since  the  medium  was  2000  m  wide.  Well,  I  now  trust  this  aspect  of  the  program.  Thenc 
work  is  to  get  to  the  Beam  program. 
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Appendix  C 

Source  Modeling  Results 


SSP  Source  Model  Sensitivity  Analysis 


In  the  process  of  studying  the  absorbing  layer  problem,  discussed  in  the  last 
section,  I  found  that  a  detailed  analysis  of  the  source  models  used  in  SSP  was  necessary. 
It  seemed  that  solutions  converged  only  after  the  grid  sizes  smaller  than  three  points  per 
wavelength.  Since  the  efforts  would  yield  preliminary  source  sensitivity  data  this  was  very 
important.  There  are  three  types  of  source  models  used  in  SSP 

Gaussian  Model 


real(i/)  =  , 


IK 

/ 

f] 

-exp 

^  kliz+z,  )^Y 

V  2 

exp 
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^  ) 

^  h 

Green's  Function  Model 


real(u(z)) 


■if 


a(r-z,)exp 
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-a{z+z,)exp 


f  iA 


where 


3.0512 

/ 

a(z)  =  (1.4467  -  0.420Uo'z') 
Normal  Mode  Model 


3.0512 


2^^ 


real(w(z))  = 


f  sin(A:,z,)  sin(V)^,, 


While  the  rule  of  thumb  in  the  under  water  acoustics  community  is  three  points  per 
wavelength  for  grid  sizes,  this  choice  is  seemingly  too  sparse  for  some  of  the  source 
models  in  SSP.  Figures  1-3  (found  at  the  end  of  the  discussion)  show  the  real  part  of  the 
solution  generated  by  the  three  source  models  for  three  frequency  cases  of  f  =  5,  10  and 
50  MHz.  The  medium  is  assumed  to  be  250  m  high,  and  the  source  is  at  100  m.  In  Figure 
1(a),  where  the  source  frequency  is  5  MHz,  and  dz  is  25  m,  the  Gaussian  startup  (rep 
circle)  is  non  zero  at  only  3  points.  The  Green's  startup  is  non  zero  at  7  points.  The 
normal  mode  is  very  different  with  its  side  lobes,  but  there  still  is  a  rough  sampling  of  the 
oscillations.  When  decreasing  dz  to  10  m,  as  in  Fig  1(b)  already  there  is  a  better  sampling. 

Figures  1(c)  and  1(d)  show  over  sampled  cues.  The  same  is  true  for  higher  frequencies  as 
show  for  f  =  10  MHz  (Figure  2,)  and  f  =  50  MHz  (Figure  3).  Since  the  oscillations 
become  so  wild  at  50  MHz  blowups  are  show  zeroing  in  on  the  interval  [75,125]  km.  ^  I 


However,  it  is  apparent  from  examine  normal  mode  startup  fields  that  for  f  =  5 
MHz,  when  dz  =  10  m  (  6  points  per  wavelength)  the  trig  functions  are  sampled  well.  For 
f  =  10  MHz,  dz  =  5  m  gives  a  similar  sampling  as  in  the  f=  5  MHz  case.  Finally  for  f  = 
50,  as  expected  the  sampling  rate  is  the  same  when  dz  =  Im.  As  for  the  Green's  function 
and  Gaussian  function  startup  fields,  one  can  calculate  the  dz  by  fixing  the  number  of 
points  needed  for  a  fixed  fall-off  rate.  That  is  if  we  let 

z  -  z,=idz 

then  the  main  contributing  exponential  is  at  1/e  when  the  argument 


k^dz^i^ 

— - =  1  for  Green's  startup 

3.0512 


and 


k^,dzY 

4 


=  1  for  Gaussian  startup. 


So  for  example  to  sample  the  main  beam  of  a  Gaussian  startup  with  5  points,  2  on  each 
side,  set  i  =  2  and  the  equation  gives 


For  f  =  5  MHz,  this  equation  gives  dz  =  9.6,  for  f  ^  10  MHz,  this  equation  gives 
dz  =  4.8,  and  for  f  =  50  MHz,  this  equation  gives  dz  ^  .96.  Checking  the  figures  (lb), 
(2b),  and  (3c),  we  see  it  is  title.  The  main  lobe  falls  to  l/e  of  its  maximum  value  by  the 
second  notch  point  on  each  side.  And  since 


Co  A 


dz  = 


X  X 
— « — , 
2;r  6’ 


CX 


this  means  that 


or  approximately  six  points  per  wavelength.  However,  from  my  test  runs  I  don't  belk'e 
this  is  a  fine  enough  mesh  for  proper  source  sampling.  I  suggest  a  minimum  of  12  points 
per  wavelength  for  convergent  solutions  (as  far  as  the  source  model  is  concerned). 

For  now  this  can  be  handled  by  tolerating  longer  run  times,  however,  for  intensive 
computation  this  might  be  too  stringent.  An  adaptive  routine  might  be  desirable.  That  is 
one  that  starts  out  with  a  mesh  for  proper  source  sampling,  then  after  a  few  steps  cuts 
back  on  the  mesh  size.  This  is  not  a  trivial  matter,  and  it  has  been  decided  to  table  that 
effort  and  study  the  real  problem  encountered. 

As  for  a  comparison  of  source  models,  it  appears  that  in  the  far  field  removed  from 
the  source,  it  really  does  not  matter  which  particular  source  one  chooses.  In  Figure  4  I 
have  plotted  the  losses  for  the  three  different  source  models.  In  the  near  field  at  ranges 
less  than  10  km  (  Fig  4(b)  and  4(c)),  the  solutions  do  disagree,  but  by  the  far-field  at 
ranges  of  20  to  25  km  (Fig  4(d)),  all  three  solutions  agree  quite  excellently.  I  thought  that 
the  absorbing  layer  might  be  influencing  the  results  so  1  ran  the  problem  again  without  it. 
These  results  are  shown  in  Figure  5.  By  ranges  of  5  km,  all  three  solutions  already  agree 
well.  (Fig  5(c)  and  5(d).).  In  Figure  6  the  frequency  has  been  raised  to  10  MHz.  Similar 
results  are  seen  with  exceUent  agreement.  Finally  in  Fig  7,  frequency  is  raised  to  25  MHz, 
and  there  is  still  excellent  agreement  between  the  solutions  in  the  far-field. 
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